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The Thermal Stability of Mass-Loaded Flows 
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Abstract. We present a linear stability analysis of a flow undergoing conductively-driven mass-loading from embed- 
ded clouds. We find that mass-loading damps isobaric and isentropic perturbations, and in this regard is similar 
to the effect of thermal conduction, but is much more pronounced where many embedded clumps exist. The 
stabilizing influence of mass-loading is wavelength independent against isobaric (condensing) perturbations, but 
wavelength dependent against isentropic (wave-like) perturbations. We derive equations for the degree of mass- 
loading needed to stabilize such perturbations. We have also made ID numerical simulations of a mass-loaded 
radiative shock and demonstrated the damping of the overstability when mass-loading is rapid enough. 
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1. Introduction 

The thermal instability of radiative media was examined 
by Field ( 1965), who considered both the presence and ab- 
sence of thermal conduction, and derived the growth rates 
of isobaric and isentropic perturbations. The first numeri- 
cal calculations of catastrophic cooling in shock heated gas 
were performed by Falle 1|19751 I1981[l . Radiative shocks 
were shown to exhibit a global overstability by Langer et 
al. Q1981[) . and have since been extensively examined (e.g., 
Chevalier & Imamura 119821 Imamura, Wolff & Duiscn 
IT351 Gaetz, Edgar & Chevalier HUM Blondin & Cioffi 
IT9531 Strickland & Blondin 113551 Walder & Folini Ir5rjfj|> . 

The interaction of radiative flows with cold embed- 
ded clouds is known to significantly modify these flows 
(see, e.g., Pittard, Dyson & Hartauist l2l)l)ll and references 
therein). However, there has yet to be an investigation 
into how mass-loading may affect the stability properties 
of such flows. This is the aim of this paper. 

In Sec. |21 we perform a linear stability analysis of 
a static medium, in thermal equilibrium, undergoing 
conductively-driven mass- loading. In Sec.[3]we present nu- 
merical models of mass-loaded radiative shocks to examine 
the suppression of thermal instability in them. We finish 
in Sec. 01 with a discussion on the conditions necessary for 
mass-loading to suppress the thermal instability in a typi- 
cal planetary nebula, and provide numerical estimates for 
the Helix nebula. 
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2. Instability in a uniform medium, including the 
effects of conductively-driven mass-loading 

The dynamics are governed by the standard hydrody- 
namic equations for plane parallel flow, 
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vP = -Cp, (3) 



together with an equation of state, 

P= -pT. 
A* 



(4) 



We have assumed that the clouds are all at their equi- 
librium radius so that if T < Tq they grow and that if 
T > T they evaporate (McKee & Cowie rHT77|) . q has the 
dimensions of a mass-loading rate per unit volume. 

The equilibrium state is characterized by p = po , T = 
Tq, v = 0, and £(pq,Tq) = 0. Assuming perturbations of 
the form 



a(z, t) = ai exp (nt + ikz), 



(5) 
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we find the linearized equations for the perturbations to 
be 



The growth of the isobaric instability (which Field 
refers to as a condensation mode) requires 



npi + p ikv! = <7 ATi/Tb, 

npoVx + ikPi = 0, 
n 7 
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and 
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Since a m is positive by definition, mass-loading always 
acts to reduce the growth rate of this instability mode 
(as Field found for conduction). However, unlike the 
corresponding equation including conduction, Eq. 1181 is 
independent of k, so mass- loading stabilizes all wave- 
lengths equally effectively against isobaric perturbations. 
Rewriting Eq.^J we find that the instability is suppressed 
if 



where C p = (d£/dp)r and Ct = (8C/dT) p are evalu- 
ated for the equilibrium state. We are left with 4 variables 
(pi, Pi, Ti, and V\) in 4 equations. The resulting disper- 
sion relation is 
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This is roughly equivalent to requiring that the mass- 
loading timescale be less than the cooling timescale of the 
medium. 

The growth of the isentropic instability requires 
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, Mass-loading is again always stabilizing, but is not equally 

The adiabatic speed of sound is c = (7P0/P0) 1 , and we effective ftt aU wavelengths in SU p pre ssmg isentropic per- 



have introduced the wavenumbers 
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turbations. Rewriting Eq. [^H] we find that the critical 
wavenumber above which perturbations are stabilized is 
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and that the isentropic instability is suppressed if 
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By introducing the non-dimensional variables 



n 

kc' 



T' 

(j'rp = UT 



Over the temperature range 5 x 10 5 <T<5x 10 7 K, 
a good approximation is C = pKT^ 1 / 2 (e.g., Kahn0976). 
This leads to greatly simplified versions of Eas.ll9land l22l 
which respectively become 
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we can write the dispersion equation in the form 

y 3 + y 2 (r' T + y(l + a p a m ) + (a T + a m j - er p )/7 = 0. (15) 

The coefficient in the y term can be removed by the in- 
troduction of the variable 



and 
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The dispersion relation then becomes 
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and to the relationship k p /kx = o p f ot — 2. 
3. Hydrodynamical calculations 
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The nature of the overstability of radiative shocks is 
known to depend on the temperature dependence of the 

0. (17) 

local cooling rate. For a power-law dependence (e.g., 
Cp oc T a ), the system is overstable for values of a below 



This equation is now in the same form as Eq. 18 in Field some critical value, a cr . Systems with a > a cr are stable. 



(1965). 



Previous numerical work has shown that a CI « 0.4 (e.g., 
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Imamura et al. 119841 Strickland & Blondin 1995), in good 
agreement with the linear stability analysis of Chevalier 
& Imamura (|1982(1 , 

To investigate the effect of mass-loading on the stabil- 
ity of a radiative shock, we have computed ID numerical 
calculations for a Mach 10 shock with the condition that 
a = 0.0. Our computations were initialized in a similar 
fashion to that presented by Strickland & Blondin (1995 ) 
for the case of an isolated planar shock, and were per- 
formed using the same hydrodynamical code, VH-1 (see 
Blondin et al. lTM]|) . 

Briefly, we assumed a steady-state radiative shock in 
the centre of the grid, with pre-shock flow from the left 
grid boundary, and cold, dense, post-shock material ex- 
iting the right grid boundary. The most common down- 
stream condition for similar calculations in the litera- 
ture is that of a reflecting wall, but a continuous outflow 
has several benefits. For instance, feedback between the 
cold dense post-shock layer and the cooling gas is prop- 
erly treated, and the problem more closely resembles the 
common occurence of an isolated shock in the interstellar 
medium. We also repeat the cautionary note of Strickland 
& Blondin (1995) that the boundary conditions play an 
important role in determining the stability of the shock. 
Therefore, it is desirable to place the shock well away 
from grid boundaries. We refer the reader to Strickland 
& Blondin (1995) for a fuller description of the code and 
initial conditions. In our runs we assumed 7 = 5/3. 

In the left panel of Fig. ^ we show the development 
of the overstability for a Mach 10 flow with a = 0.0 and 
qo = 0.0 (i.e. no mass-loading). The pre-shock density and 
velocity are set to unity. At the beginning of the simula- 
tion, corresponding to the top of the plot, the flow is ini- 
tialized as closely as possible to the steady state solution 
for a radiative shock with no mass-loading. The oversta- 
bility is excited from weak perturbations produced by the 
numerical mapping of the steady state solution onto the 
computational grid. Linear growth was observed for about 
the first 40% of the run, after which the amplitude of the 
oscillations saturate and the system enters the nonlinear 
regime. The behaviour shown in the left panel of Fig. [T] is 
in excellent agreement with Fig. 3 of Strickland & Blondin 
( 1995) where the overstability for a Mach 40 flow is dis- 
played. 

In the right panel of Fig. ^ we show the evolution 
of a mass-loaded radiative shock with qo = 150.0 and 
A = 5/2. The process of mass- loading increases the density 
and thermal pressure, and decreases the temperature and 
velocity of the flow. The pressure increase causes the shock 
position to initially expand away from the dense cooling 
layer. However, the enhanced density leads to more rapid 
cooling, and thermal pressure support is soon lost. This 
causes the shock to reverse its direction of motion and 
to fall towards the cold dense layer. The shocked gas is 
then repressurized and this cycle is repeated. With each 
cycle, the amplitude of the oscillations decreases as the 
mass-loading damps this instability. 



With a = 0.0, £ = 1, /.t/R = 1 and the post-shock 
values p = 4, T = 0.193, Eg. 1191 shows that we require 
qo > 13 for mass-loading to begin to damp the overstabil- 
ity. In practice, we find that we need values appreciably 
larger than this for strong damping, although the oscil- 
lation amplitude does decrease slightly for values of qo 
near 13. There are a couple of reasons why we need higher 
values of qo than suggested by Eq. ^5] The principal rea- 
son is due to the fact that adding mass to the shocked 
gas decreases the cooling timescale while simultaneously 
increasing the mass-loading timescale. Hence, while the 
initial addition of mass to the immediate post-shock flow 
may be at a rate such that the mass-loading timescale is 
less than the cooling timescale, the difference between the 
two timescales will diminish with time and may even be 
reversed. 

Secondly, the large dynamic range in temperature 
which exists for high Mach number shocks (e.g., for a 
Mach number of 40, the immediate post-shock temper- 
ature is almost 3 orders of magnitude greater than the 
ambient temperature) leads to a smaller value of T avg /T ps 
than is the case for lower Mach number shocks. Here T avg 
is the average temperature over the cooling length of the 
shock, and T ps is the immediate post-shock temperature. 
We find that as we reduce the shock Mach number, the 
value of qo needed for significant damping more closely 
matches that from Eq. EH 

While we have not been concerned with the ioniza- 
tion state of the gas in this model, it is worth noting that 
the radiative cooling rate can be greatly enhanced when 
neutral gas is introduced into hot plasma, and that the 
ionization energy can also be significant under such con- 
ditions (Slavin et al. lTM3|) . 

4. Discussion 

As a relevant example, the above results can be applied 
to planetary nebulae (PNe). Such objects often display 
clumps which appear to mass-load the shocked wind of 
the central star. The clumps were presumably formed by 
instabilities in the atmospheres of the red supergiant stars, 
prior to the evolution of the central star into a hot object 
with a fast wind. In the following we derive relationships 
for the number density of clouds and their total mass, with 
the condition that mass injection from the clouds is able 
to suppress the thermal instability. 

Since the central stars of PNe appear to fall into 3 
groups (e.g., Kudritzki et al. 1997), with winds spanning a 
large range of mass- loss rate (M ~ 10 -9 — 10 -6 M Q yr _1 ), 
we choose not to focus on an individual nebula. Instead, 
we note that the typical thermal pressure in PNe, wind- 
blown-bubbles, and starburst superwinds is comparable 
to the pressure in clumps embedded within them, P/k ~ 
10 7 cm -3 K. In PNe, this must balance the ram pres- 
sure from the wind of the central star. The post-shock 
number density is then n = 4 x 10~"' 'v^ 2 'k / /i(P '/ 'k)j = 
55vf 2 (P/fc) 7 cm , where V7 is the wind speed in units 
of 1000 km s" 1 and (P/k) 7 is P/k in units of 10 7 cm" 3 K. 
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g cm 3 s 1 , where 
1.33 x 10- 19 /fi 2 and a r- 1 



In the last step, and throughout the remainder of these 
calculations, we set n = 10~ 24 g. The post-shock temper- 
ature, T — 1.36 x 10 5 v 2 K. 

To suppress the thermal instability by mass-loading, 
Eq. d gives q > 10~ 32 ?;f 7 (P/fc)f. 
we have adopted A — - , 

(cf. Kahn 11575)1 . and substituted for T and n. The mass 
evaporation rate from a single clump is rh — 2.75 x 
10 4 T 5 / 2 i? pc (30/lnA) gg -i (Cowie&McKeelMB, where 

T is the temperature of the ambient surroundings, and R pc 
is the radius of the clump in parsecs. Since InA ~ 30, rh ~ 
1.9 x 10 17 Wyi? pc gs _1 , and the required number density 
of clumps is n c = q /m > 1.6 x 10 6 uf 12 R-^(P/k)^ pc~ 3 . 

Observations of the Helix nebula (NGC 7293) indi- 
cate that the density within the clumps, N c ~ 10 6 cm -3 , 
while they are unresolved on spatial scales of ~ 10 -3 pc 
(O'Dell & Handron UMB) . Adopting R pc = 10~ 3 , v 7 = 2, 
and (P/fc)7 = 1 yields a number density of clumps of 
n c > 4 x 10 5 pc -3 . Since the radius of the Helix nebula is 
~ 0.2 pc, we then expect ~ 1600 clumps, with a total mass 
of ~ 0.1 M©. While these estimates are in good agree- 
ment with observational inferences from the Helix nebula, 
small changes in v-j can greatly affect these values. For 
the evaporation to be smooth on a spatial scale of l/fc m , 
we require that the clumps are smaller than L = 2n/k m , 
i.e. L < 1.5 X 10~ 2 Vj(P/k)i 1 pc. With the above parame- 
ters, L < 1 pc, and is consistent with the derived number 
density of clumps. 

Acknowledgements. We would like to thank the referee, John 
Raymond, for timely and constructive comments. JMP would 
also like to thank PPARC for the funding of a PDRA posi- 
tion. This research has made use of NASA's Astrophysics Data 
System Abstract Service. 



References 

Blondin, J.M., Cioffi, D.F., 1989, ApJ, 345, 853 

Blondin, J.M., Kallman, T.R., Fryxell, B.A., Taam, R.E., 1990, 

ApJ, 356, 591 
Chevalier, R.A., Imamura, J.N., 1982, ApJ, 261, 543 
Cowie, L.L., McKee, C.F., 1977, ApJ, 211, 135 
Falle, S.A.E.G., 1975, MNRAS, 172, 55 
Falle, S.A.E.G., 1981, MNRAS, 195, 1011 
Field, G.B., 1965, ApJ, 142, 531 

Gaetz, T.J., Edgar, R.J., Chevalier, R.A., 1988, ApJ, 329, 927 
Imamura, J.N., Wolff, M.T., Durisen, R.H., 1984, ApJ, 276, 
667 

Kahn, F.D. 1976, A&A, 50, 145 

Kudritzki, R.P, Mendez, R.H., Puis, J., McCarthy, J.K., 1997, 
in IAU Symp. 180, Planetary Nebulae, eds. H.J. Habing 
and H.J.G.L.M. Lamers (Dordrecht: Kluwer), 64 
Langer, S., Chanmugam, G., Shaviv, G, 1981, ApJ, 245, L23 
McKee, C.F., Cowie, L.L., 1977, ApJ, 215, 213 
O'Dell, C.R., Handron, K.D., 1996, AJ, 111, 1630 
Pittard, J.M., Dyson, J.E., Hartquist, T.W., 2001, A&A, 367, 
1000 

Slavin, J.D., Shull, J.M., Begelman, M.C., 1993, ApJ, 407, 83 
Strickland, R., Blondin, J.M., 1995, ApJ, 449, 727 
Walder, R., Folini, D., 1996, A&A, 315, 265 



time 




Fig. 1. Spacetime df|@i\^4i^(lc45il)9 : tiY e snoc k for M = 
10, a = 0.0 (cf. Bia 3 in Staaklaitti ^ ^^^ [TM^i . The 

tick marks on the vertical axis are in mnits of 25L c /vq 
where L c is the cooling length of the shock (see Strickland 
& Blondin [1995). The width of the horizontal axis is 
1.88L C . The supersonic flow enters the grid from the left, 
and the cooled postshock gas flows off the grid to the 
right. The grayscale shows the density of the gas, with 
lighter shades signifying diffuse gas and darker regions de- 
noting higher density. In the left panel, <7o — 0, & n d the 
system exhibits linear growth for the first ~ 40% of the 
run. In the right panel, where go = 150, the initial con- 
ditions are far from equilibrium for a mass-loaded shock, 
and the system exhibits powerful oscillations, though these 
are quickly damped as expected from the linear stability 
analysis in Sec. [21 



